Spinful bosons in an optical lattice 
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We analyze the behavior of cold spin-1 particles with antiferromagnetic interactions in a one- 
dimensional optical lattice using density matrix renormalization group calculations. Correlation 
functions and the dimerization are shown and we also present results for the energy gap between 
ground state and the spin excited states. We confirm the anticipated phase diagram, with Mott- 
insulating regions of alternating dimerized 5=1 chains for odd particle density versus on-site 
singlets for even density. We find no evidence for any additional ordered phases in the physically 
accessible region, however for sufficiently large spin interaction, on-site singlet pairs dominate lead- 
ing, for odd density, to a breakdown of the Mott insulator or, for even density, a real-space singlet 
superfiuid. 

PACS numbers: 32.80.Pj, 71.35.Lk, 03.75.Mn, 75.40.Mg 



I. INTRODUCTION 

Strongly interacting systems which are at the forefront 
of study in theoretical condensed matter physics can be 
realized with cold atomic gases in optical lattices. Spinor 
atoms in optical lattices constitute a novel realization of 
quantum magnetic systems. These new realizations have 
several advantages compared to their condensed matter 
counterparts. One is precise knowledge of the underlying 
microscopic model, another is the possibility to control 
the parameters of the lattice Hamiltonian, a third is the 
absence of impurities in the system. Due to the small- 
ness of the scattering length compared to inter-particle 
separation, a gas of degenerate alkali-metal atoms is con- 
sidered as a weakly interacting gas Ij . However this is no 
longer true if the atomic scattering length is changed by 
means of a Feshbach resonanceQ or if an optical lattice 
created by standing wave laser beams is used to confine 
the atoms in the minima of the lattice potential, which 
strongly enhances the effects of the interactions. 

For the case of the optical lattice the Mott insulator- 
superfluid quantum phase transition was first demon- 
strated in the seminal experiment of Greiner et al.|3|. for 
®^Rb. This experiment has led to an enormous research 
activity in the field of cold atoms confined in optical lat- 
tices. In most experiments on cold atoms a magnetic field 
is used to create a trap that confines the atoms in the 
lattice. This means that the spins of the atoms are fully 
polarized so that the atoms behave as spinless particles. 
Recently a new experimental setup that uses an optical, 
instead of magnetic, trap has been developed0,0]. Using 
this type of confinement atoms with different spin po- 
larizations are trapped and the scattering of the atoms 
become spin dependent 0, 0, 0- The spin interaction 
may be either antiferromagnetic or ferromagnetic in its 
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structure depending on the scattering length which is 
material specific. For sodium ^■^Na, the interaction is 
antiferromagnetic 6] . 

A model for spinful bosons in optical lattices has been 
studied before by a number of groups "t*. 's', fiol HTl IT^ 
■It^j . and a phase diagram has been predicted. For integer 
number of atoms per lattice site there are two regimes. 
One, where the kinetic energy is dominating over the 
potential energy, has a superfiuid groundstate and the 
other, where the interaction energy is the most impor- 
tant, has a Mott insulating groundstate. For non-integer 
particle density the superfiuid groundstate is always pre- 
vailing. For the one-dimensional case the predicted phase 
diagram has a spin-dimerized phase in the insulating re- 
gions with odd particle density and on-site spin singlets 
in the insulating regions with even density. This phase di- 
agram was recently confirmed numerically However, 
for the study of the order parameter the authors' resorted 
to a mapping of the boson model to a spin modelf^. 
By doing so the dimension of the local Hilbert space is 
reduced from 20 or more (depending on the maximum 
number of bosons per site) to three, and accordingly the 
calculations are less time and memory consuming. How- 
ever, this mapping is only applicable to the insulating 
regions with odd density, and is valid in the limit of very 
strong atom-atom repulsion. 

Our aim with this article is to study the spinful model 
directly. The correlation functions and excitation ener- 
gies are presented for the first three insulating regions. 
We also compare results for the insulating systems with 
those for the superfiuid phase. All calculations in this ar- 
ticle are done using density matrix renormalization group 
(DMRG) with open boundary conditions. 

The outline of the paper is as follows; In Sec. |n] the 
model and the mapping from the spinful Hamiltonian to 
the spin-one chain is presented. In Sec. IIIII some details 
about the calculations are given. In Sec. II VI we present 
particle-particle and spin-spin correlation functions for 
the system. In Sec. the dimerization in the spin-spin 



correlation is obtained and from this it is concluded that 
the third Mott lobe is dimerized and that the first lobe 
most probably is dimerized. In Sec. IVII the energy gap 
to excited spin states is investigated. It is apparent that 
the second Mott lobe has on-site singlets and that the 
gap decays with a universal behavior when the superfluid 
phase is approached. We also show that the odd density 
insulating systems have a small spin gap and that this 
energy gap is non-zero for all spin-interaction strengths. 
Section IV ill contains a discussion on the conditions un- 
der which bound on-site singlet phases occur. This gives 
several new states, including a singlet insulator and a su- 
perfluid phase of condensed singlet pairs (SSC) d, 13 • bi 
this latter phase the tunneling of individual atoms is sup- 
pressed and only singlet pairs tunnel^G^]. The last section, 
Sec. I Villi is devoted to a summary and conclusions. 



II. THE MODEL 

In most experiments on optical lattices a magnetic trap 
is used to confine the atoms. Due to the spin dependent 
interaction between the atoms and the trap only one of 
the spin species is confined in the lattice. This s^tem is 
very well described by the Bose-Hubbard model[l^]. The 
phase diagram for the spinless Bose-Hubbard model was 
discussed in the seminal work by Fisher et al.jlTj using 
a scaling theory and renormalization-group calculations. 
They also derived the exact phase diagram within mean- 
field theory, i.e. for an infinite-range hopping model. The 
phase diagram for the one-dimensional model has been 
obtained numerically with a high precision by Kiihner et 
al. J^] using the DMRG method. The phase diagram 
has Mott insulating regions (Mott lobes) and superfluid 
regions. A sketch of the phase diagram is shown in the 
left panel of Fig. 

When an optical trap is used instead of a magnetic 
trap, atoms with all spin polarization are trapped, with 
alkali atoms having hyperflne spin F — 1. The scattering 
between two atoms takes place in the spin-zero and the 
spin- two sectors. Since the scattering length is different 
in these two sectors the on-site repulsion becomes spin 
dependent, resulting in the following Hamiltonianj^ Q, 

(1) 

+ ^(i7n,(n, -1) + JS^) . 

i 

Here bf^ and bi^a are the creation and annihilation oper- 
ators for bosons with z-component of spin a — 0, ±1, 
in the lowest Bloch band localized on site «, = 
Scr ^ta^i,c is the local density, t is the hopping inte- 
gral for atomic wavefunctions between different lattice 
sites, U is the usual on-site Coulomb repulsion. Si — 
Tlia a'^tJ^a a'^i+i a' s^ui Operator for spin one 

particles at lattice site i, with ^' being the usual spin- 
1 matrices. J is the spin dependent interaction whose 
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FIG. 1: A sketch of the phase diagram of the Bose-Hubbard 
model f]3 . The left panel shows the Mott lobes for the spinless 
Bose-Hubbard model and the right panel shows the spinful 
phase diagram {J/U = 0.1). 

value and sign depend on the material. For ^^Na, J is 
positive, i.e., the spin interaction is antiferromagnetic, 
and J/U « 0.040. In this paper we focus mainly on the 
parameter range appropriate for ^^Na. 

The insulating regions with an odd density form an ef- 
fective spin-one chain, as the spins on each lattice site will 
combine to the lowest possible spin. This maps the in- 
sulating odd density region for the spinful Bose-Hubbard 
model onto the bilinear-biquadratic spin-one chain with 
HamfltonianQ, S E 

H = ^(cos6l(Si • Sj) + sin6l(Si • Sj f) (2) 

(y) 

where tan(0) = jiz^jjiji £ [— 37r/4, — 7r/2], for the first 
Mott lobe. Similar expressions should exist for higher 
(odd) Mott lobes. 

The mapping is obtained by making a lowest order, 
i.e. second order, perturbation expansion in t of the 
Hamiltonian Eq. (Q. To second order only pairwise in- 
teractions between atoms on neighboring sites are gener- 
ated. For the first Mott lobe, i.e., in the insulating state 
with exactly one boson per site, the nearest-neighbor in- 
teractions are always of the form given by Eq. (|2Jl. How- 
ever, the derivation of the 6 dependence on J and U is 
only valid in the limit t U. Away from the limit t <^U 
but still in the insulating phase, higher order terms have 
to be included which will, in addition to renormalizing 
the nearest-neighbor couplings, add spin couplings be- 
yond nearest-neighbor. 
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Despite many years of study, the phase diagram of the 
spin 1 chain is not fully resolved. For 6 = 0, the an- 
tiferromagnetic Heisenberg model, the spectrum has a 
finite Haldane gapJT^, whereas at 9 = ±7r the system 
is ferromagnetic. Actually, both phases have a finite ex- 
tension in parameter space. Ferromagnetism exists for 
— TT < 9 < —3tt/A and tt/2 < 9 < tt while the massive 
Haldane phase exists for ~n/4 < 9 < tt/A. At the lower 
end of this interval, — 7r/4, the gap vanishes, but opens 
again for 9 < — 7r/4 and a massive dimerized state with 
spontaneously broken translation symmetry is found. For 
the special point 9 = — 7r/2, the ground-state energy, gap, 
correlation length and the dimer order parameter can be 
calculated exactly [23. The unresolved issue is whether 
the dimerized phase extends all the way to 9 = —3tt/A 
or if there exists another phase in between the dimerized 
and the ferromagnetic phase near 9 = — 37r/4. 

Chubukov|2H. studying fluctuation effects near the 
end point of the ferromagnetic phase, concluded that 
the dimerized and the ferromagnetic phase are sepa- 
rated by a disordered phase, a gapped non-dimerized 
nematic. According to Chubukov a direct transition be- 
tween the dimerized and ferromagnetic phases is very un- 
likely since completely different symmetries are broken in 
these phases. Chubukov claimed that the dimer order pa- 
rameter and the gap vanish simultaneously at a 9c above, 
but close to, 9 = — 37r/4. For 9 < 9c the gap opens up 
and closes again at 9 — — 37r/4 whereas the dimer or- 
der parameter is zero in this interval. Buchta et al. re- 
cently performed highly accurate DMRG calculations [2^ 
for this region. Their results indicate that the dimer 
phase prevails down to 9 = — 37r/4, although they can- 
not rule out that a non-dimerized phase exists in an ex- 
tremely small interval of 9 close to 6* = — 37r/4. A recent 
preprint by Lauchli, Schmid and Trebstl23!|, use a strong- 
coupling series expansion which gives a vanishing gap for 
-37r/4 < 9 < -0.6777, but they also present DMRG re- 
sults that are consistent with those of Buchta et al.[2^ 
showing a small but finite gap in this region. Porras 
et al.|24| introduced a new numerical algorithm and ap- 
plied it to the same problem. They found indications 
of a quantum phase characterized by nematic quasi-long- 
range order in the interval — 37r/4 < 9 < — O.Ttt. However 
the accuracy of this calculation is difficult to quantify; 
their results are also consistent with either the correlation 
length being longer than the size of the chains considered, 
or a gap that is smaller than the numerical accuracv|24j. 
Also, Rizzi et al.0| concluded from DMRG calculations 
that there is no intermediate nematic phase, although 
they found indications that a tendency towards nematic 
ordering is enhanced as the dimer order parameter goes 
to zero. 

In the odd density insulating state with more than one 
boson per site there is an additional constraint in order 
to arrive at the Hamiltonian Eq. configurations with 
spin on individual sites higher than one have to be ne- 
glected. Matrix elements for scattering into such states 
are of the order {nt)'^/U and their energy is set by J, 
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FIG. 2: (Color online) The population of the single-site ba- 
sis states for an average density of one, two, A'^ = 2L -\- 2 
and three particles. The probability for the different sectors 
is presented as a function of the on-site repulsion U/t, for 
J /U = 0.05, in the first three panels. The last panel shows the 
third Mott lobe as a function of J /U , with U/t = 15. The dif- 
ferent sectors, (n, S), presented in the figure are: star(0,0), di- 
amond(l,l), triangle up(2,0), triangle down(2,2), square(3,l), 
circle(3,3), triangle right(4,0). 



where n is the number of bosons per site. Therefore in 
order to obtain the bilinear-biquadratic Hamiltonian, the 
condition is that nt < (t/J)^/^ 7]. 

The insulating regions with an even density behave 
completely differently compared to those with odd den- 
sity. For an even number of spins per site, the particles 
form bound states of singlet pairs resulting in a rather 
large energy gap to the excited states. Note however 
that unless the spin interaction is very large, the even 
density superfluid shows no such real-space pairing. 

The spin interaction affects the critical interaction 
strength for the transition from the Mott insulating to 
the superfluid phase. In the Mott lobes with even parti- 
cle densities, i.e. the insulating regions with on-site spin 
singlets, the insulating regions are stabilized by the spin 
interaction. On the other hand the Mott lobes with an 
odd density are weakened by the spin interaction, the 
density fluctuation is enhanced since the particle hop- 
ping is the process that carries the spin-interaction |25l| . 
In the right panel of Fig. a sketch of the new phase 
diagram is shown together with the phase diagram for 
the spinless model (left panel). 
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III. METHOD 

The method we have used is the Density Matrix Renor- 
mahzation Group[2^|23|) adapted to exactly handle the 
non-Abelian SU(2) symmetry |2g. The starting point of 
this formulation is the Wigner-Eckart theorem 29] , which 
states that the matrix elements of SU (2)-invariant oper- 
ators factorize into a component that is independent of 
the geometry (i.e., the quantization axis), and a geomet- 
rical factor (Clebsch-Gordan coefficient). In particular, 
we have the relationship 

(/m'|Al3ljm) = Ci^C X (j'||AW|b-) , (3) 

where (j'||A['^1||j) is the reduced matrix element of the 
rank k tensor AI*^!. From this relation, all of the oper- 
ations required for the DMRG (tensor products, wave- 
function transformations etc.), can be written in terms 
of the reduced matrix elements only, with an appropriate 
6j or 9j coupling coefficient 28, 29]. 

The use of symmetries to improve the efficiency 
of DMRG calculations was recognized from the 
beginning 2G| , and presumably arose from the Wilson's 
use of spin in NRG|30t. the historical forerunner of 
DMRG. The utility of Abelian symmetries (typically 
U{1), such as particle number, or z-component of spin), 
is that it imposes a sparseness constraint on the opera- 
tor matrix elements; given a basis labeled by the parti- 
cle number, for example, the matrix elements of some 
irreducible operator {n'\A^-^^\n) are non-zero only for 
n' = n + N. By storing only the matrix elements that 
are permitted by symmetry to be non-zero, the memory 
use and computational cost is greatly improved. Note 
however that a priori this has no effect on the accuracy, 
and in principle for a given basis size the obtained matrix 
elements and wavefunction will be identical whether the 
symmetry constraint is used or not. This is in contrast 
to the case of non-Abelian symmetries, where Eq. © 
implies a reduction in the basis size needed to represent 
the Hamiltonian matrix elements. That is, a single basis 
state of total spin j requires 2j + 1 states to represent in 
a t7(l) basis. This results in a proportional reduction in 
the basis size m required for a given accuracy. Depend- 
ing on the spin of the block basis states this reduction 
factor typically ranges from around 3 in the Mott insu- 
lating phase where the block spins tend to be minimal, to 
8 or so in the superfluid phase. Since the computational 
cost of the algorithm is proportional to m^, this results 
in orders of magnitude improvement in the efhciency. 

For these calculations, we typically used 350 basis 
states (equivalent to approximately 1000 — 2500 basis 
states of a U{1) calculation). For the finite size scaling 
analyses we used up to 300 lattice sites for the phases 
with total number of bosons A'' equal to the number of 
lattice sites L, up to 120 sites for N = 3L and up to 
70 sites for the N = 2L superfluid. We fixed the maxi- 
mum number of bosons per site to be 5, except for the 
N = L superfluid where we found four bosons per site to 



be sufficient, and the N = L insulator where we found 3 
bosons per site to be sufficient. This is justified by the 
occupation number shown in Fig. [21 where the average 
number of bosons in the different single site basis states 
are shown. 



IV. CORRELATION FUNCTIONS 

In this section we discuss the calculated hopping cor- 
relation functions and the spin-spin correlation func- 
tions. In the superfluid the hopping correlation func- 
tion r(j) = {bf^bi+j_a) decays with a power-law as we 
expect p^. The low-energy energy physics of the super- 
fluid phase is described by a Luttinger liquid and the 
decay of the hopping correlation function Tij) ^ j^^^"^ 
where K is the Luttinger liquid parameter |3l|. Since 
we use open boundary conditions, measures have to be 
taken to reduce the effects of the boundary. The most 
important effect is oscillations in the local density. In the 
superfluid phase they show the characteristic Luttinger 
liquid power-law decay away from the edge of the sys- 
tem. These density fluctuations will affect the hopping 
correlation function r(j). Since we are interested in the 
properties of the infinite system we reduce the effects of 
these fluctuations by averaging over pairs of {bfa^i+j.a) 
for fixed j. We discard contribution from the ten lattice 
sites closest to the boundary where the Friedel oscilla- 
tions are largest. Beyond this distance the numerical 
behavior of F is no longer sensitive to the boundary os- 
cillations. The right panel of Fig. O shows results for 
the hopping correlation function for two different integer 
densities displaying superfluid behavior. 

In the Mott insulating phase the hopping correlation 
function r(j) is decaying exponentially 17], where the de- 
cay rate is determined by how close we are to the phase 
boundary of the Mott phase. As one would expect, the 
boundary effects decay exponentially in this region. The 
left panel of Fig. O shows results for the hopping correla- 
tion function for three different integer densities. 

In Fig. 21 the spin-spin correlation functions for differ- 
ent systems are presented for Mott insulating and su- 
perfluid regions. The spin interaction strength is set to 
J/U — 0.05. The figure demonstrates that the correla- 
tion function is rather independent of the phase for the 
odd density systems. It decays in the same way for all 
the odd density systems and also for the superfluid sys- 
tem with an even density. The insulating system with a 
density of three atoms per site is dimerized, for this in- 
teraction strength, which is the reason for the oscillating 
structure of the spin-spin correlation. 

For the insulating system with an even density the 
spin-spin correlation is exponentially decaying, consistent 
with the picture of on-site singlets. 
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FIG. 3: The (6^6)-correlation function. The left hand panel 
is in the insulating regime {U/t — 10, U/t = 10, U/t = 15 for 
N = L, N = 2L, N = 3L respectively), the right hand panel 
shows the superfluid {U/t = 1 for both N — L and — 2L), 
log- log plot (Inset: same data with a log- linear scale). We fix 
J/U = 0.05 in all cases. 




FIG. 4: The spin-spin correlation for N = L, both insulating 
and superfluid, N = 2L both insulating and superfluid and 
N = 3L insulating, J/U — 0.05. The clear oscillation in the 
correlation function for the third Mott lobe arises from the 
large dimerization. 

V. DIMERIZATION 

One of the questions that we try to resolve in this arti- 
cle is the nature of the spin order in the first Mott lobe. 
The phase is, as mentioned above in Sec. ^ believed 
to be dimerized. The standard way to characterize the 
dimerization is by using a dimer order parameter which 
usually is defined as the difference in bond energy in the 
middle of the chain 

D = - (i/j-i,,). (4) 

In the spinful Bose-Hubbard model this dimer order pa- 
rameter is given by {b1jbj^i—b^_.^hj) since all other terms 



FIG. 5: The absolute value of the difference in neighboring 
correlation functions in the flrst Mott lobe {U/t = 10, J/t = 
0.5) for different system sizes. The dimerization is calculated 
up to the middle of the chain. 



in the Hamiltonian are on-site. This difference is very dif- 
ficult to measure since it is very small, see Fig.O The fig- 
ure shows the absolute value of the mentioned difference 
as a function of the lattice site, j, for a Mott insulating 
system. In the figure j runs between 1 and the middle of 
the chain. The dimerization is calculated for a couple of 
different system sizes to illustrate the boundary effects. 
The bond energy is indeed oscillating between different 
lattice sites but the values are so small that it is difficult 
to tell whether or not the oscillation will be non-zero in 
the thermodynamic limit. The open boundary conditions 
break the translational symmetry which is necessary to 
see a dimerized order parameter. For periodic boundary 
conditions the groundstate would be a linear combination 
of the two possible degenerate dimerized configurations 
and translational symmetry would be restored for any 
finite L. 

The dimerization in the spin-spin correlation functions 
for the Mott insulating system has also been studied, the 
result is presented in Fig. |31 In this figure the dimeriza- 
tion in the spin-spin correlation is shown as a function of 
the lattice site, for a few system sizes. 

(Color online) The dimerization in the spin-spin inter- 
action for different interaction strengths in the first Mott 
lobe is shown in Fig. El In the figure the spin-spin dimer- 
ization at the middle of the lattice is shown as a function 
of the system size. There is a dimerization when the 
spin interaction becomes large but it is hard to draw any 
conclusions about the thermodynamic limit for small in- 
teraction strengths. We have done a finite size scaling of 
the points with the following fitting fimctionp^l 

D = Da + dN-^ exp(- (5) 

The infinite size values for the first Mott lobe are pre- 
sented in Fig. [7| The dimerization goes to zero alge- 
braically with an exponent ^ 6, i.e., Dq cx {J/U)^. This 
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FIG. 6: (Color online) The dimerization in the spin-spin cor- 
relation function as a function of the inverse lattice size for 
the first Mott lobe {U/t = 10), for different strength of the 
spin interaction. The solid lines are the finite size scaling. 
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FIG. 7: The dimerization in the thermodynamic limit as a 
function of spin interaction strength. The dimerization in the 
first Mott lobe decays algebraically with a fitted exponent 
~ 5.7 and there is no suggestion that it will be non zero 
before J/U = 0. In the third Mott lobe the dimerization is 
much larger. 



FIG. 8: (Color online) The difference in energy between states 
with S = 2 and the ground state, S = 0, as a function of the 
inverse system size for U/t — 10 and different spin interaction 
strengths. The values for the second Mott lobe (denoted '-I-') 
are scaled by 1/100. 

nite size scaling is larger. As was mentioned in Sec. Hll 
also for the third Mott lobe the spinful Bose-Hubbard 
model could be approximately mapped to the bilinear- 
biquadratic spin chain. However in this case the con- 
straints for the perturbation expansion in t to be valid 
are more severe, configurations with spin on individual 
sites higher than one have to be neglected leading to the 
condition nt <C (J7J)^/^. For the dimerization in the 
third Mott lobe presented in Fig. [3 this condition reads 
J/U 3> 0.04. As configuations with spin higher than one 
come into play we see a possible cross over in the behav- 
ior and a down-tm'n in the curve for Dq for J < 0.01. 
However still we find no indications of Dq being zero for 
J larger than zero. The third Mott lobe is definitely 
dimerized for J/U = 0.005 and we conclude that there 
is no extra phase between the dimerized phase and the 
ferromagnetic in the third Mott lobe. 



VI. SPIN GAPS 



is in agreement with a calculation on the bilinear bi- 
quadratic chain' 14]. The curve shows no indication of 
a non-zero dimerization before J = 0. The values in the 
figure are obtained for U/t = 10 and for U/t = 3. The 
fact that the curves ior U/t = 3 and U/t = 10 fall on top 
of each other shows that the dimerization is present in 
the entire Mott lobe. This conclusion cannot be deduced 
from the bilinear-biquadratic chain, since the mapping is 
only valid for small hopping. 

In Fig. the dimerization in the spin-spin correlation 
for the third Mott lobe as a function of spin interaction 
strength is also shown. The values are infinite size values 
obtained in the same manner as for the first Mott lobe. 
However, the system sizes are smaller so the error in fi- 



In this section, the difference in energy between the 
ground state, which is a spin singlet, and the first excited 
spin state, which has spin 5 = 2, is studied. This differ- 
ence is calculated for several system sizes and is shown 
as a function of the inverse system size in Fig. |S| The 
energy gap to excited states in the first Mott lobe has 
been calculated before with DMRG using the mapping 
to the bilinear-biquadratic spin chain 14J and also in the 
limit of zero hopping 32] . 

From the figure it is obvious that the decay of the en- 
ergy gap, with respect to the inverse system size, for the 
Mott insulating state with two particles per site behaves 
differently to that of the other cases. The gap is large 
which is consistent with the on-site-singlet picture and 
approaches the thermodynamic limit as L~^|32|. 
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FIG. 9: The spin gap in the thermodynamic limit for the 
second Mott lobe. In the left panel the values are shown as 
a function U/t. In the right panel the lines are shifted by a 
constant Au/t along the x-axis to show the similarity in the 
decay of the gap for the different interaction strengths. 
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FIG. 10: (Color online) The difference in energy between 
states with 5 = 2 and the ground state, S = 0, as a func- 
tion of the inverse system size for U/t = 10 and different spin 
interaction strengths the solid lines are finite size scalings. 



We have studied how the energy gap in the thermo- 
dynamic hmit for the second Mott lobe approaches the 
superfluid value as the on-site repulsion, U, is decreased. 
The result is shown in Fig.El for J/U = 0.05, J/U = 0.1 
and J/U ~ 0.15. To obtain the values in the figure cal- 
culations are done for different interaction strengths and 
system sizes. The infinite size values are found by fitting 
second-order polynomials to the data points. In the fig- 
ure the infinite size values are presented as a function of 
the on-site repulsion for three different strengths of the 
spin-spin interaction. In the right panel of Fig.|^the ab- 
scissa (U/t) of two of the curves in the left panel has been 
shifted by a constant (A^/f), different for each of them, 
to show that the curves can be collapsed on top of each 
other. From this it is clear that the energy gap decays 
in the same manner for the three interaction strengths 
when the superfluid phase is approached. 

We have further studied how the gap approaches zero 
for the odd density insulating regions when J/U ap- 
proaches zero. In Fig. ^1 we show for the first Mott 
lobe aX U/t = 10, the gap as a function of system size, 
for a variety of different spin couplings. The symbols 
mark the numerical data and solid lines show the extrap- 
olation to the thermodynamic limit using a second order 
polynomial fit in 1/L. 

In the first Mott lobe the gap in the thermodynamic 
limit is very small and decays to zero as the interaction 
approaches zero, see Fig. 1111 In the figure, results for 
U/t = 10, i.e. far inside the Mott lobe, are shown to- 
gether with values for U/t — 3 which is fairly close to the 
phase transition. Results for the third Mott lobe and the 
superfluid phase are also shown in the figure. There is 
no indication of that the energy gap becomes zero before 
J = 0. The four curves in Fig.^jcan best be fitted with 
a power law decay and nothing suggests an intermedi- 
ate phase in between the dimerized phase and the J < 
ferromagnetic phase. 
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FIG. 11: (Color online) The energy gap between the ground 
state and the first spin excited state in the thermodynamic 
limit, as a function of the spin interaction for odd density 
systems. The system with U/t = 1 is superfluid. 



In Fig. ^] we show the energy per site as a function of 
the magnetization. We show examples for both insulat- 
ing and superfluid systems. The spin excitation spectra 
for the first and the third Mott lobe and the superfluid 
regions follow the parabolic form typical of Heisenberg 
spin systems [3^, 



E{S)=EQ + kS{S+l). 



(6) 



The behavior of the spin excitations for the second 
Mott lobe is rather different, see Fig. ^1 In this case, 
the energy increases linearly with S and with a steep 
slope, a result which is consistent with the picture of the 
system being composed of on-site singlets. The energy 
to break another singlet is independent of the number of 
already broken singlets. 

The energy gap to excited spin states measures the 
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FIG. 12: The energy density as a function of the magneti- 
zation M = S/ L. The open circles for the second Mott lobe 
have been scaled by a factor 20 to fit on the figure. The fit for 
these points is a straight line. All other fits are of the form 
E(, + kS{S+l), J/U = 0.05. 



magnetization energy and as expected from the on-site 
singlet picture the second Mott lobe has a higher mag- 
netization energy. 



VII. SINGLET CONDENSATES 

In this section we comment on some of the proposed 
phases for the spinful Hubbard model. Zhou and Snoek 
have suggested that a phase exists in between the su- 
perfluid and the Mott insulating phases for even parti- 
cle density [ll- In this hypothesized phase, the particles 
would form a condensate of tightly bound on-site singlet 
pairs. However, we see no evidence for this phase for 
any physically realistic magnitude of spin interaction. In 
Fig. 13 we show the occupation of each site basis state, 
for J/U ~ 0.05. For this value of the spin interaction, it 
is clear that there are no values of the Coulomb repul- 
sion U where only even number of bosons is favored. This 
also follows from a simple energy argument, that suggests 
such a phase is unlikely: Since the spin interaction energy 
is smaller than the charging energy it is costly to increase 
the density on a site by two. The energy gain by keeping 
the spins coupled to a singlet is much smaller than the 
energy penalty for the on-site Coulomb repulsion. 

For large coupling the system behaves similarly to a 
spinless Bose-Hubbard model populated with spin sin- 
glets. This limit is tested and indeed for J /t — 10 and 
N = L the system is essentially composed oi L/2 empty 
sites and L/2 double occupied sites, see Fig. ^| The 
limit J ^ U,t could localize the singlets to give a crystal 
phase:6j, but we have not investigated this scenario. In 
this figure we also show, as a function of U, the popula- 
tion in the different on-site number and spin sectors for 
N — 2L particles per site for the spin interaction strength 
J/t =10, and also the superfluid system with N = 2L + 2 



1 



a) 



"I □ 



0.01^^ A 

<5 O 



0.0 



a iii C ] 



b) 



6 o 



01 0.1 u/t 

c) 



O.lr 



0.01^ 



A 
O 



A 

A 



1 U/t 10 



d) 



0.1 



0.01 



O.lr 



V V 37 



O A 



U/t 



- <> o o o ^ 

I A 

0.001 ^ t , t , , 



0-1 u/t 



FIG. 13: The probability to have a certain number of parti- 
cles and spin on a lattice site is presented for different systems 
and different density and spin (n,S)sectors, circle(0,0), dia- 
mond(l,l), square(2,0), triangle up(3,l), triangle down(4,0). 
a.)N — 2L and J/t = 10: once U and t are comparable in size a 
singlet condensate forms. b)A'^ = 2L and J/U = 1: for large I] 
this is an insulating spin-singlet state, seen from the compar- 
atively high probability in the n = 4 sector. c)A'^ = 2L-I-2 and 
J/U = 10: spin singlet condensate. d)Af = L and J/t = 10: 
in this spin singlet condensate the lattice is half-filled with 
spin singlets. 



for J/U — 10. As J/U is reduced, the charging energy 
becomes equal to the spin energy and the singlets are no 
longer bound. 



VIII. CONCLUSION 

We have studied the ID spinful Bose-Hubbard model 
using DMRG. By examining the energy and correlation 
functions we have resolved the behavior of the model for 
all of the parameter ranges that are physically relevant 
for optical trapping experiments. 

Much evidence is presented showing that the insulating 
system with an even density is well described by the on- 
site spin-singlet picture. The spin gap in the thermody- 
namic limit is rather large, paired with a short correlation 
length for the spin-spin correlation function. The energy 
increases linearly with the spin of the system, showing 
that the energy required to break another singlet is given 
by the energy gap and is independent of the number of 
already broken singlets, thus the singlets are strongly lo- 
calized. 

It has been more difficult to determine the nature of 
the other insulating regions. From the particle-particle 
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correlation function it is concluded that there is a super- 
fluid and an insulating phase just as for the spinless case. 
It has previously been suggested that these insulating re- 
gions are dimerized. The dimer order parameter which 
uses differences in bond energy is difficult to measure in 
numerical simulations since it is rather small. Instead 
the oscillations in the nearest neighbor spin-spin corre- 
lation were used to show that it is likely that the odd 
density per site Mott lobes are always dimerized. This is 
a relatively strong result for the third Mott lobe, where 
the dimerization order parameter is quite large. For the 
first Mott lobe our results are consistent with the S — 1 
Heisenberg calculations of Buchta et al. 22] and Rizzi et 
al.^3|- In general the numerical accuracy of calculations 
in this regime will be much better for the Heisenberg 
model compared with the spinful Bose-Hubbard model, 
so our result does not contribute substantially to the ar- 
gument over the absence (or otherwise) of the Chubukov 
nematic phase. However what is clear is that the dimer- 



ization in the first Mott lobe shows scaling behavior with 
an exponent ~ 6 throughout the entire Mott lobe. We 
found no evidence for novel superfluid phases; in partic- 
ular the on-site paired superfluid suggested by Zhou and 
Snock 8] is absent in this model. 

The energy gaps to spin excited states were studied, 
showing that the odd density (superfluid and insulating) 
and even density superfluid phase follow the form of a 
parabolic spectrum, proportional to S{S + 1). 
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